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We use a lattice gas cellular automata model in the presence of random dynamic scattering sites 
and quenched disorder in the two-phase immiscible model with the aim of producing an interface 
dynamics similar to that observed in Hele-Shaw cells. The dynamics of the interface is studied as 
one fluid displaces the other in a clean lattice and in a lattice with quenched disorder. For the 
clean system, if the fluid with a lower viscosity displaces the other, we show that the model exhibits 
the Saffman-Taylor instability phenomenon, whose features are in very good agreement with those 
observed in real (viscous) fluids. In the system with quenched disorder, we obtain estimates for 
the growth and roughening exponents of the interface width in two cases: viscosity-matched fluids 
and the case of unstable interface. The first case is shown to be in the same universality class of 
the random deposition model with surface relaxation. Moreover, while the early-time dynamics 
of the interface behaves similarly, viscous hngers develop in the second case with the subsequent 
production of bubbles in the context of a complex dynamics. We also identify the Hurst exponent 
of the subdiffusive fractional Brownian motion associated with the interface, from which we derive 
its fractal dimension and the universality classes related to a percolation process. 

PACS numbers: 47.56.+r, 47.11.Qr, 68.05.-n 


I. INTRODUCTION 

The dynamics of Newtonian fluids is governed by the 
Navier-Stokes equations [T]. In fact, a number of dis¬ 
tinct flow behavior has been identified, such as laminar 
or turbulent, depending on boundary conditions and the 
Reynolds number, Re = VL/v, where V (L) is some 
characteristic velocity (length) and ir is the kinematic 
viscosity of the fluid. A special case of wide interest is 
that of fluid dynamics in porous media mis]. 

In a microscopically disordered media the flow is very 
complex and its macroscopic behavior can be homo¬ 
geneous or heterogeneous [2]: in the former case, the 
system displays size-independent transport properties, 
while in the latter the system is better described by a 
position-dependent permeability Further, at very 

low Reynolds number with the characteristic length es¬ 
timated from, e. g., the average size of the solid ob¬ 
stacles, a proper volume average of the porous media 
strongly indicates, in agreement with experimental ob¬ 
servation, that the single-phase flow is well described by 
Darcy’s law Hi. In addition, if dispersion (diffusion) 
effects are relevant, an advection-dispersion (convection- 
diffusion) equation must be included in the description 
of a variety of phenomena, such as fluid flow and solute 
transport [5], and miscible viscous fingering Hi. On 
the other hand, as Re increases, violations of Darcy’s 
law have been reported due to inertial effects, described 
by the Forchheimer equation |ini lllj , and very interest¬ 
ing non-Newtonian behavior, with power-law permeabil¬ 
ity and data collapse in a broad range of Reynolds con¬ 
ditions [12]. Very recently, the stability analysis of two- 
phase buoyancy-driven flow, in the presence of a capillary 
transition zone, was investigated in detail m- 


A large body of work on fluid dynamics in porous me¬ 
dia clearly indicates that several powerful contin¬ 

uum and discrete approaches have been put forward to 
treat many cases of interest. In the continuum version 
both analytical and direct numerical integration of the 
pertinent dynamic equations have proved very efficient 
HHH [i EHUI, while several concepts from percola¬ 
tion, growth models, random walk and fractal geome¬ 
try have been very useful in dealing with discrete mod¬ 
els H [m US], which are alternative approaches heav¬ 
ily based on numerical simulations and supplemented by 
continuum stochastic equations m- Among the discrete 
models widely used in the literature, we shall emphasize 
the lattice gas cellular automata. 

Lattice gas cellular automata (LGCA) is a family of 
computational models whose particles have velocities de¬ 
fined by a discrete set and collide with each other on sites 
of an hexagonal lattice; the resulting velocity configu¬ 
ration obeys mass and momentum conservation dlHIS]. 
Making use of the Chapman-Enskog expansion, it was 
demonstrated that in two [HHini and three [20] dimen¬ 
sions (2D and 3D) the average values of the pertinent 
variables in a macroscopic scale satisfies the incompress¬ 
ible Navier-Stokes equations; the 3D case requires mi¬ 
nor additional considerations. More recently, 2D-LGCA 
flows on curved surfaces with dynamical geometry were 
also proposed |21] . with motivation in describing a vari¬ 
ety of phenomena. 

Ghallenging aspects in the description of flow in porous 
media are associated with the complex geometry of the 
pore space and the region near the fluid-fluid interface in 
the multiphase case. The LGGA model handles the no¬ 
slip boundary condition in the fluid-solid interface in a 
very simple manner through the local bounce-back rule. 
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which is easily implemented even for the complex ge¬ 
ometries arising from the pore space |22j . In particular, 
the LGCA single phase flow in 2D porous media sat¬ 
isfies Darcy’s [53] and Forchheimer ’s m laws, and the 
permeability of real systems can be estimated mn EH 
from images of their microgeometry. We also remark 
that, in a clean channel, Poiseuille flow is observed [53] . 
Two-phase flows in 2D are also modeled through the 
LGCA by the use of colored particles and an additional 
collision-color conservation law |25|. This model leads 
to a well-defined interface between the two fluids, with 
a finite surface tension, thereby allowing the observation 
of Laplace’s law and the phenomenon of phase separa¬ 
tion [25]. Beyond that, it was shown [26] that this model 
captures the expected interface hydrodynamic fluctua¬ 
tions associated with capillary waves, and exhibits dy¬ 
namic scaling exponents |16j in agreement with those of 
the Kardar-Parisi-Zhang equation [221 EH], apart from 
logarithmic corrections for the largest systems studied 
[26j . The macroscopic properties of two-phase immisci¬ 
ble fluids in 3D porous media were described [29] using 
lattices built from rock’s microgeometry. Further, the 
flow of two-phase immiscible fluids in the presence of an 
amphiphile species in 2D m and 3D m porous media 
were also considered. 

We also stress that much attention has been devoted to 
the development of the lattice-Boltzmann method (LBM) 
[351133] . This variant of the LGGA method [13] can be 
simplified through a linearized collision operator [331 EH] , 
both in 2D and 3D. The LBM was also used in the study 
of single- [361 El] and two-phase [38ll45] flows in artifi¬ 
cial and real 3D porous media; as well as in the study 
of single-phase flow in multiscale porous media [3]- Fur¬ 
ther, the LBM was adapted [33] to model single-phase 
flow in a 2D porous medium, in which case the action of 
the solid space is introduced as linear (Darcy’s law) and 
non-linear (Forchheimer’s law) effective body-force terms 
in the generalized Navier-Stokes equations [44]. On the 
other hand, viscous finger phenomenon in 2D channels 
was studied through the LBM for two miscible fluids 
[45j and two immiscible fluids with zero |45j and finite 
[46] surface tension. We also mention that the LBM was 
adapted to describe the hydrodynamics of 3D liquid crys¬ 
tals m, and very recently used [35] to study the flow of 
a nematic liquid crystal in microfluidic channels. 

In the following we digress on the statistical properties 
and the dynamics of an interface separating two immis¬ 
cible fluids flowing in a Hele-Shaw cell, or in 2D and 3D 
systems modeling these cells, which is the main subject 
of this work. The dynamics of growing interfaces has 
stimulated intensive theoretical and experimental inves¬ 
tigations motivated by its own challenging fundamental 
aspects and potential technological applications as well 
[IS] [IS]. In this context, the displacement of a fluid by 
another during the flow in a porous medium has been 
considered as a representative example in the study of 
the dynamics of roughening interfaces [THIITS]. The main 
experimental tool in such studies is the Hele-Shaw cell 


Si ED], which mimics the porous media and have un¬ 
veiled a rich variety of phenomena. The cell is formed by 
two plates separated by a little gap of size b. The veloc¬ 
ity of a viscous single-phase flow through the gap, aver¬ 
aged through the direction perpendicular to the plates, 
satisfies Darcy’s law with permeability 6^/12. For an im¬ 
miscible two-phase flow in a horizontal cell, if the less 
viscous fluid displaces the more viscous one, the inter¬ 
face is unstable, and a linear stability analysis predict the 
occurrence of one-finger or multifinger patterns |49L 150] . 
More recently, very detailed experimental aspects of this 
phenomenon [3IHS3], and theoretical studies based on 
vortex-in-cell method [S3], random-walk techniques and 
methods from complex analysis [S3], have been put for¬ 
ward; which also include phase-field model [55] , potential 
flow analysis [56] of radial fingering [57] [S5] and rigorous 
arguments [59]. 

In the context of flow in clean Hele-Shaw cells, a sec¬ 
ond variant of the LGCA method was formulated [60| [BT] 
shortly after the original proposal [13 [n], in which case 
the Darcy’s law obtains without the need of the detailed 
microscopic geometry of the medium. In fact, through 
the introduction of random scattering sites, in a stochas¬ 
tic context, the viscous dissipation of the single-phase 
flow in a Hele-Shaw cell was investigated using this mod¬ 
ified LGCA method. On the other hand, in the pro¬ 
posed model for immiscible fluids [65], the authors used 
dynamic random scattering sites (a time-dependent ap¬ 
proach) and find that the interface fluctuations exhibit 
typical random-walk behavior, with a roughness expo¬ 
nent 1/2; however, the fluctuations are not associated 
with hydrodynamic capillary waves, since the kinetic 
rules [55], and other assumptions, differ from those of 
Ref. [55]. Notwithstanding, Saffman-Taylor fingering in¬ 
stability [491 EQ] is observed through this modeling [53] : 
the interface is unstable either if the fluids differ in pres¬ 
sure or in the density of dynamic random scattering sites. 
On the other hand, we mention that a two-phase model 
similar to that of Ref. [25], but with only local update 
rules [53] and with the inclusion of stochastic scattering 
sites, was suggested [55] to describe multifractal prop¬ 
erties of fingering patterns associated with the Saffman- 
Taylor instability in a clean Hele-Shaw cell. The viscosity 
contrast was introduced through the use of distinct col¬ 
lision rules for each fluid. The authors find that for high 
surface tension the Hausdorff (fractal) dimension is in the 
interval (1.41 - 1.46), while for moderate surface tension 
the fractal dimension is found to be 1.75 (see discussion 
below). 

We now discuss highlights of viscous-fingering phenom¬ 
ena in Hele-Shaw cells with quenched disorder. In the 
case of a circular cell, fluid invasion at the center of the 
cell exhibits fractal fingering structure for high capillary 
number Ca = [iVja, where fj, = pv is the dynamic vis¬ 
cosity, p is the density of the fluid, and cr is the inter¬ 
face tension [55]. On theoretical grounds, for simplicity, 
the viscosity of the invading fluid is considered negligible 
(infinity viscosity ratio), thus allowing a map [53155] of 
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the growing front dynamics onto the diffusion-limited ag¬ 
gregation (DLA) process [69], characterized by a fractal 
dimension Df « 1.71. Notwithstanding, the measured 
[66j fractal dimension of the fingering structure for high 
Ca, Df = (1.62 - 1.64)±0.04, is very close to the frac¬ 
tal dimension of the percolation cluster backbone [TH] . 
Df = (1.61 - 1.62) ± 0.02; it also satisfies the inequal¬ 
ity, 1.3 < D^JIj^ < Dyl < 1.7, suggested [71] to be 
valid above the percolation threshold in the context of 
the analogy between DLA, percolation and viscous fin¬ 
gering. However, at low Ca the fractal structure is de¬ 
scribed by invasion percolation with trapping [Z21E3], 
Df « 1.82, which should be compared with the fractal 
dimension of ordinary m or invasion [ZHESKin] per¬ 
colation at the threshold: Df — 91/48 « 1.89. We also 
remark that the diffusion front of random walkers on 2D 
lattices is a fractal object of dimension Df = 1.75, in 
close connection with the percolation cluster hull EH; 
in 3D lattices the front is a dense object with D = 3 
EH]. Further, in rectangular cells, two-case studies have 
been reported m- (i) water (wetting fluid) displacing oil 
(nonwetting fluid); (ii) oil (nonwetting fluid) displacing 
water and glycerol (wetting fluid), with viscosity ratio 
M = ~ 200 in the two cases. The authors find 

that in the first case (imbibition process, with the finger 
width A ^ Ca“ 2 ) ^ the finger width is much larger than 
the pore size (viscous fingering), whereas in the second 
case (drainage process), the finger width has the same 
order of the pore size (capillary fingering). In this con¬ 
text, a dynamical transition has been identified |80| . On 
the other hand, self-affine (anisotropic) fractal interfaces 
m I81j from immiscible displacement in porous media 
were observed when the displacing fluid is more viscous 
and more effectively wets the medium (e. g., water or 
glycerol displacing air in a Hele-Shaw cell filled with glass 
beads). In these experiments the scale-dependent rough¬ 
ness, \V{L) = AL°', where W{L) is the interface width, 
a « 0.73 is the roughness exponent, and A ~ Ca “2 
[82j . while large variations with time and flow condi¬ 
tions suggest a more complex dependence on Ca [55] . 
with the effective value of a varying over a wide range 
(0.65 - 0.91). In the latter experimental study, connec¬ 
tion with the dynamics of random-field Ising model was 
proposed, although in this case a crossover from a = 0.8 
at small scales to a = | at larger scales was found [84] , in 
agreement with experiments on wetting invasion |85L I86j . 
We remark that, in the context of fractal growth in hy¬ 
drodynamic dispersion through random porous media, 
direct integration of the advection-diffusion and linear 
stokes equations shows that the following steady flows 
are obtained [57] according to the value of the Peclet 
number, Pe = where Dm is the diffusion constant: 
(i) for Pe = 00 , the fractal morphology of the spread¬ 
ing dye (tracer) is characterized by a fractal dimension 
close to Ddla] (ii) for Pe = 0, a self-affine fractal inter¬ 
face is identified with scaling consistent with a growth 
process as predicted by the KPZ equation [27]; (hi) for 


finite Pe, the behavior is characterized by the competi¬ 
tion between advection and diffusion effects. Lastly, a 
series of very detailed experimental studies in Hele-Shaw 
cells with disorder have investigated the interface scaling 
and fractal properties under combined viscous, gravity, 
and capillary effects [88ll93j . In particular, for immisci¬ 
ble displacement of viscosity-matched fluids in 2D porous 
media [51], the fractal structures of the front are char¬ 
acterized by using the box counting algorithm, which 
implies Df = 1.33 ± 0.05, consistent with the fractal 
dimension of the external perimeter in invasion perco¬ 
lation [M], and the density-density correlation function 
(data collapse for several values of Ca), which implies 
Df = 2 — a = 1.42 ± 0.05 and /3 = 0.8 ± 0.3, where /3 is 
the growth exponent {Wl ^ t^)- In addition, the authors 
find that the fractal dimension of the structure left be¬ 
hind the front is Df = 1.85, consistent with invasion per¬ 
colation with trapping. Further, a detailed analysis of the 
two-phase flow of immiscible fluids in disordered porous 
media allowed the authors [H] to propose a quite unified 
view of the macroscopic transport properties. Indeed, the 
characteristic crossover scales between fractal regimes, i. 
e., from capillary to viscous fingering, are discussed us¬ 
ing box counting, where the box side size I ranges from 
the cell width Ly (cell extent L^) down to the pixel size. 
They thus find: D — 1.00 for I > Ly, Df = 1.60 for 
Ly > I > a/Ca, and Df — 1.83 for I < a/Ca, where a is 
the pore size. We stress that the two mentioned fractal 
dimensions are consistent with those of 2D percolation 
cluster backbone and invasion percolation with trapping, 
respectively. 

The main goal of our work is to study the interface dy¬ 
namics between two immiscible fluids in a Hele-Shaw cell 
geometry in clean systems and in systems with quenched 
disorder; both cases are investigated through the LGCA 
model |25] in the presence of dynamic random scatter¬ 
ing sites [63] . In fact, our approach combines the best 
features of the LGCA model, namely, its microscopic 
kinetic rules [25j . interface hydrodynamic fluctuations 
effects [55], and dynamic random scattering sites [55] . 
thereby allowing a nice description of the Saffman-Taylor 
instability in a clean Hele-Shaw cell, in very good agree¬ 
ment with experimental observations [IH] [SD] [HI [S3] . In 
addition, the inclusion of quenched disorder allows us to 
consider the competition between viscous and capillary 
effects at the pore level, thus causing the roughening in¬ 
terface to exhibit very interesting self-affine features. 

The work is organized as follows, in Section II we 
briefly review the LGCA models from which we build the 
two phase LGCA with scatterers. In Sec. HI we present 
our results for the interface dynamics in a clean system, 
exhibiting the Saffman-Taylor instability, and in a sys¬ 
tem with quenched disorder, in which case we also study 
scaling behavior, critical exponents, and fractal dimen¬ 
sion of the interface. Finally, a discussion on the main 
results and concluding remarks are presented in Section 
IV. 
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II. LGCA MODELS 
A. Single Phase LGCA 

The basic model reviewed in this section is named the 
random-collision 7-velocity [12]. Time and space are dis¬ 
crete and the particles can lie only on the sites of a regu¬ 
lar triangular lattice. Moreover, its velocities can assume 
only seven possible values, one is null (particle at rest) 
and six have magnitude of one lattice unit (unit of length) 
per time step (unit of time) in the direction of the lat¬ 
tice edges, Ci = (cos(i7r/6),sin(i7r/6)) with i = 1,...,6, 
where a lattice unit is the distance between two sites. 
Not more than one particle per site in one direction is 
allowed. In each time step the outcome of particle colli¬ 
sions are such that total momentum and particle number 
are conserved at each site of the lattice, after each colli¬ 
sion the particles hop to a neighboring site in accord with 
its velocity. In a collision, the model takes into account 
all possible randomly changes in the velocity allowed by 
the conservation laws; the collision is local, i.e., depends 
only on particles at the same site. 

In this model, it is natural to define the following hy¬ 
drodynamic variables: 

6 

(1) 

i^O 

the density of particles per site at position x and time t; 

6 

( 2 ) 

i=l 

the momentum density per site, where Ni is the prob¬ 
ability of a particle at position x and time t to have a 
velocity in the direction i, considering Nq in Eqs. 0 
and ([^ as the probability of a particle to be at rest. 
These probabilities are obtained by making space-time 
averages. Densities per unit area, p and g, are given by 

p = p/(V3/2) and g = g/(V3/2) (3) 


marked as solid sites. In our simulations the force is 
introduced by increasing the momentum at a randomly 
chosen site, if such a change is possible, until the total 
desired momentum increase is reached in each time step. 

B. Immiscible Two Phase LGCA 

The Immiscible Lattice Gas (ILG) was introduced by 
Rothman and Keller |25j and mimics the dynamics of two 
immiscible fluids in two dimensions. This model is sim¬ 
ilar to what we have described above except that there 
are two kinds of particles: red and blue. As before, we 
allow not more than one particle with the same veloc¬ 
ity at the same site, even if they have different colors. 
In the propagation, the particle maintains its identity 
and in a collision, in addition to the conservation of mo¬ 
mentum and particle nnmber, the number of particles of 
each color is also conserved at each site. Nonetheless, 
not all collisions are allowed; in fact, they depend on the 
neighborhood and in the way that particles of the same 
color tend to cluster. More specifically, the site config¬ 
uration could be represented by 14 Boolean variables, 
{ro,ho,... ,rQ,be) each one representing the presence of a 
red or blue particle in the direction i (again, the direction 
0 is associated with a particle at rest). The effect of the 
neighborhood on the result of a collision at the site x is 
introduced through the definition of the vector quantity 
color gradient: 

F(x) = J2ciY^[rj{x + Ci) - bj{yi + Ci)]; (6) 

* j 

the vector color flux: 

q[r(x), b(x)]) [r,(x) - 6,(x)], (7) 

i 

and by considering that the resultant site configuration 
maximizes the quantity: 

q(r',b')-F, (8) 


In systems not far from equilibrium, with smooth vari¬ 
ations in space and time, Ghapman-Enskog expansions 
allows one to obtain the following dynamic equation [m- 

dl: 

^+ 5 (/o)(u-V)u=--Vp-f h, (4) 

ot p p 

where u is the average flow velocity, v is the kinematic 
viscosity, f is the body force per unit area and 


9{p) 


14 7 - 2p 
M 7-p ■ 


(5) 


Boundary conditions are easily implemented [TTI - 
m and no-slip boundary condition is introduced by 
bouncing-back the particles when they arrive at sites 


and conserves the number of each kind of particle and the 
total momentum. In the above equation, prime indicates 
the configuration after a collision. If there are more than 
one configuration that maximizes Eq. Q, one of them is 
randomly chosen. 

Through these rules, an initially mixed system will sep¬ 
arate into regions with only red or blue particles, depend¬ 
ing upon the density of particles and the relative concen¬ 
tration of color [191 [25] • In general, if the density is above 
2 particles per site and the relative concentration is not 
too high or too low, the separation occurs. In a phase 
separated system it is verihed by simulations that there 
is a discontinuity in pressure which obeys Laplace’s law 

[BUS]: 


Pl-P2= CFK 


(9) 
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FIG. 1. (Color online) Profile of the x-component of the mo¬ 
mentum per site, g^, for the single phase LGCA with scat- 
terers in a flow with the same boundary condition as the 
Poiseuille flow (indicated in the fignre). The density of scat- 
terers, Na, is indicated, the density per site is p = 3.5 and the 
body force per site is / = 1/(256 x 64) « 6.1 x 10“®. Solid 
lines are fittings of the data to Eq. ( |12[ ). Inset: dissipative 
factor Ua as a function of Na from the presented profiles. 


where pi and p 2 are the pressure in regions with different 
colors, cr is the surface tension and k is the curvature of 
the interface. The value of cr varies with the density 
of particles and for the density used in this paper it is 
approximately 0.37. In regions where there is just one 
kind of particle, the dynamics is the same as that of the 7- 
random-velocity model, but in the interface the dynamics 
changes. 


C. Single Phase LGCA with Scatterers 

This model was introduced by Hayot et al. |60) with 
the goal of reproducing the fluid dynamics ruled by 
Darcy’s law. The scatterers provide momentum dissi¬ 
pation similarly to the solid region of porous media or 
the plates of Hele-Shaw cells. Here the scatterers are 
introduced in the random-collision 7-velocity model. In 
each time step a fixed number of sites is chosen randomly 
as solid scattering sites, thus inverting the velocity of the 
scattered particles and introducing a damping term in 
Eq. g. The number of scattering sites must be very 
low compared with the number of particles in the system. 
Furthermore, if the density is chosen to be 3.5 particles 
per site g(p) = 0, thereby eliminating the advective term 
in Eq. In this manner, we obtain m the following 
steady state equation 


OsU = —-Vp -I--I- -f, (10) 

P P 

where as is positive and proportional to the density of 
scatterers [SD]. For a lattice gas in the xy plane subject 
to a force density / and a constant pressure (p) gradient 


also in the x direction, Eq. (10) reads: 


asU^ = + vV‘^Ux + -f, ( 11 ) 

p ax p 

The effect of the dissipative term can be understood by 
comparing the velocity profile for the boundary condi¬ 
tions of a Poiseuille flow with the parabolic profile which 
is obtained in the absence of this term. For this bound¬ 
ary condition the fluid is confined between two lines of 
solid sites at p = 0 and y = L, which implies Ux{0) = 0 
and Ux{L) = 0, due to non-slip boundary conditions in 
the solid sites, and the following solution to Eq. 11 [60] : 


, ^ ^ fi cosh[r(y-(L/2))]) 

1' - ™h(rL/2) 1 • 

( 12 ) 

with r = 

In Fig. 12 we present the velocity profile for a simula¬ 
tion (5^ of the single phase LGCA model with scatter¬ 
ers, adopting periodic boundary conditions (PBC) in the 
X direction, for the indicated values of density of scat¬ 
terers, TVs, homogeneous density of force per site /, as 
well as data for a simulation without scatterers. Since 
the applied force is homogeneous, the density (p) and 
the pressure (p) remain constant on average (^ = 0); 
therefore, the data of the simulations are well fitted by 
Eq. (12). As we can see, the presence of scatterers flat¬ 


tens the velocity profile in comparison with the parabolic 
one, except near the walls [13 HU HT]. For Ns = 0, we 
obtain v = 0.196, which is in accordance with reported 
results [TS] for the value of v for the random-collision 
7-velocity model. Following the procedure described in 
Ref. |5D], we obtain v = 0.197, 0.176, 0.207 and 0.196 for 
Ns = 0.001,0.002,0.004 and 0.008, respectively. Also, 
we find the dissipative factor = 1.95As, which is 
near the theoretical estimate from the microscopic model 
through the use of the Boltzmann approximation m-- 
as = 2fVs, but greater than the value observed for the 
6-velocity model [50]: Ug ~ l-lNg. 

In the region where the velocity is flat, the Laplacian 
term in Eq. 0 can be neglected, and this equation 
reduces to Darcy’s law 


(13) 


k f dp 

Ux = - { + f 

p \ ax 


where p = vp \s the dynamic viscosity of the fluid and k 
is the permeability of the porous medium; for fluid flow in 
a Hele-Shaw cell, k = 6^/12, where b is the gap between 
the two plates of the cell. This implies 


k 

P 


1 

OsP’ 


(14) 


which allows a physical interpretation of the underlying 
role played by the scatterers in the model. For homoge¬ 
neous pressure, Eq. (13) can also be written in terms of 


the momentum and force densities per site: 

- - h 

9x f ; 


(15) 
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FIG. 2. (Color online) Component x of the momentum per 
site, g^, as a function of the body force per site, /, for p = 3.5 
and periodic boundary conditions in the x and y directions, 
for the values of Ns indicated. Soli d lines are fittings to a 
Darcy’s law equation form, Eq. (15l. Inset: (From Fig. 1) 
y/k = 2asp/%/3 as a function of Ns calculated from Eq. (141 
and the estimated value of Us from Fig. 1, 


In Fig. [^we present the simulation of a flow on a lattice 
without solid sites and with periodic boundary conditions 
in the x and y directions and without any obstructed 
sites. The simulations were performed using 10® time 
steps and the results for the average momentum are ob¬ 
tained discarding the first 20000 time steps. A linear 
relation between / and Qx or Ux is obtained. The values 
of ^ = a^p shown in the inset of this figure are consis¬ 
tent with those obtained from the simulation in a channel 
with as « 2Ns, in agreement with the estimate from the 
Boltzmann approximation |60] . 

Therefore, the approach with scatterers is quite re¬ 
markable since it allows the unified description of Hele- 
Shaw cell experiments with fluid dynamic viscosity tuned 
by M = {^)ctsP, or the flow in a general 2D porous me¬ 
dia, which may include quenched disorder effects, with an 
average homogeneous permeability k fixed by Eq. (14|. 

An interesting situation is the possibility of introducing 
scattering sites in the ILG model, in particular for red 
and blue fluids having distinct scattering rates, which 
can be interpreted as two fluids with different viscosities. 
This situation can be realized and shall be discussed in 
next section for an homogeneous lattice and a lattice with 
quenched disorder. 


III. IMMISCIBLE LGCA WITH DYNAMIC 
SCATTERERS 


This model is built by introducing random dynamic 
scattering sites in the ILG model. In regions of the lattice 
where there is just one kind of fluid, the flow dynamics 
is the same as the LGGA model with scatterers and the 
change occurs only on the interface. In a static situation, 



FIG. 3. (Color online) Snapshots of a simulation in which the 
red (light gray) fluid with a density of 0.001 scatterers per site 
displaces the blue fluid, which has a density of 0.008 scatterers 
per site. The forcing rate is 0.0001 and acts from left to right. 
The width of the system is 500 -\/3/2. Simulation time step 
increases from top to bottom. 



FIG. 4. (Color online) Snapshots of a simulation in which the 
red (light gray) fluid with a density of 0.001 scatterers per site 
displaces the blue fluid, which has a density of 0.008 scatterers 
per site. The forcing rate is 0.0005 and acts from left to right. 
The width of the system is 500 %/3/2. Simulation time step 
increases from (a) to (f). 


simulations as in Ref. [5S] confirm that Laplace’s law 
holds and the surface tension has a value nearly equal to 
the one found in the ILG model (« 0.37). In particular, 
if a flow is established in a system with distinct scatter¬ 
ing rates for the red and blue fluids, we can observe a 
Saffman-Taylor-like instability [49l [50] in the interface, 
as described below. 


A. Interface Instabilities 

The dynamics of the interface between two immiscible 
fluids as one displaces the other in a Hele-Shaw cell can 
be discussed by considering the average velocity of flow 
and the relation between the viscosities of the two flu- 
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FIG. 5. (Color online) Snapshots of a simulation in which the 
red (light gray) fluid with a density of 0.001 scatterers per site 
displaces the blue fluid, which has a density of 0.008 scatterers 
per site. The forcing rate is 0.0005 and acts from left to right. 
The width of the system is 2000 \/3/2. Simulation time step 
increases from (a) to (f). 


ids [IHHSS]- Three general behaviors for the interface are 
observed as these parameters are varied: (i) If a more 
viscous fluid displaces a less viscous one or a less viscous 
displaces a more viscous with a low average velocity, the 
interface is flat; (ii) if a less viscous fluid displaces a more 
viscous one with intermediate values of the average veloc¬ 
ity, a single stable finger is observed; (iii) for sufficiently 
higher values of the average velocity a multifinger dy¬ 
namics sets in. If fluid 2 displaces fluid 1 in a cell with 
width Ly the finger patterns observed are governed by 
two dimensionless parameters [5ni3T|: the viscosity con¬ 
trast 


Ml M2 
fcl k2 , 
Ml I M2 > 
ki ' k2 


and a dimensionless surface tension 


B = 


2a 1 


(16) 


(17) 


where V is the average velocity of the flow. Further, from 
a linear stability analysis of the problem, it is possible to 
show that the perturbation with the wavelength [151 HO] 


\m, — 27r L/>i 




3(7 


(t-f)^ 


(18) 


presents the maximum growth rate; however, all pertur¬ 
bations with A > Am/'s/s are unstable. In the case of a 
Hele-Shaw cell, with /i 2 << /ii, Xm ~ TrbjVC^, where 
Ca = fjX^/a is the capillary number. 

Turning our attention to the ILGS, with the red fluid 
displacing the blue fluid, we observe similar regimes to 
those discussed above and the parameters A, B and Xm 
can be written as 


— iVs,r 


(19) 


B = 


2a 1 

Lll.95{Ns,b-Ng,r)pV' 


( 20 ) 


Am = ‘2^'^Ly\ -Vb = 2tt 


3a 


1.95{Ns^b- Ns,r)pV' 


( 21 ) 


where we have used (/i/fc) = with Og = 1.95Ng, in 
Eqs. ( |I^ , ( [l7| ) and (18); moreover the subscripts r and 
b refer to the red and blue fluids, respectively. 

First, if the scattering rate of the red fluid is higher 
than that of the blue fluid, the interface is flat as in (i) 
above. If the red fluid has a lower scattering rate, we ob¬ 
serve the three situations exemplified in Figs. [^|4]and[^ 
which are obtained in a lattice under the following con¬ 
ditions: the first and last lines are made of solid sites; 
the density of both fluids is 3.5 particles per site, peri¬ 
odic boundary conditions are applied in the horizontal 
direction with the color of the particle changing when it 
crosses the boundary sites, as made in m to simulate 
a fluid invasion process. The viscosity contrast for the 
three simulations is A = 0.78. For a finite forcing rate, 
the velocity is not constant and depends on the fraction 
of each fluid; its average is calculated by considering the 
time needed to the most advanced point of the interface 
to meet the right boundary of the lattice. In particu¬ 
lar, in Fig. we show a simulation with a low forcing 
rate (=0.0001) and an average velocity V « 4.9 x 10“^ 
presenting the formation of one finger, which is our ver¬ 
sion of the so-called Saffman-Taylor instability |S3]. The 
value of Am ~ 400 implies Ly/Xm ~ 1-1, in agreement 
with the simulation. For a higher forcing rate (=0.0005), 
and an average velocity V « 2.2 x 10“^, we observe in 
Figs. 1^ and a multifinger dynamics for the interface 
between the two lattice gases, as in regime (iii) above, 
and in excellent agreement with the observed two-fluid 
dynamics in a Hele-Shaw cell [52] for the two-finger (Fig. 

and multifinger (Fig. cases. Indeed, in these two 
cases. Am ~ 190 implying Ly/Xm ~ 2.3 (« 9.1) for the 
simulation presented in Fig. El (Fig. m. Moreover, the 
agreement between our predictions for the number of fin¬ 
gers, shape and time evolution, as shown in Figs. [^|^and 
in is remarkable |5IH53]. We also emphasize that the fin¬ 
ger dynamics is quite complex; indeed, as time increases, 
some fingers can give rise to both bubble formation and 
finger amalgamation [see Figs, [^(e)-(f)]. The mentioned 
agreement is, in fact, the real reliability test of our ap¬ 
proach. 

While we do not provide a rigorous mapping between 
the immiscible two-fluid flow in a Hele-Shaw cell and the 
ILGS model, the simulations discussed above suggest the 
close relation between the two systems. 


B. Interface dynamics with quenched disorder 

The dynamics of the interface between two fluids were 
observed in a great variety of experiments in Hele-Shaw 
cells with some kind of disorder [T6| [73l l88H^ . The dis¬ 
order can be introduced by the spatial variation of the 
gap between the plates of the cell |HS] or by the introduc¬ 
tion of randomly distributed circular obstacles [551IM] ■ 















FIG. 6. (Color online) Distribution of beads radii. The data 
is obtained averaging each frequency trough 100 realization 
of initial configuration and is well fitted by a exponential dis¬ 
tribution. 


In both cases the dynamics of the interface may obey scal¬ 
ing laws, but the values of the exponents are usually dis¬ 
tinct [86l |89] . Discarding the effect of gravity, the inter¬ 
face is unstable to finger formation if the defending fluid 
is the higher viscous one, while viscous forces stabilize 
the interface in the opposite case. In the unstable situa¬ 
tion, the finger patterns can be controlled predominately 
by capillary forces (acting in the scale of the pores), or 
by viscous forces, in which case crossover length scales 
between the pole-scale structure and macroscopic scales 
can be defined [^ . 

Inspired by these experiments, we consider now the 
simulation of ILGS model for the displacement of the 
blue fluid by the red fluid in a system with a random 
distribution of circular beads (see Fig. [^. We use no- 
slip boundary conditions at solid sites, periodic boundary 
conditions in the vertical and horizontal directions, but in 
the horizontal direction the color of the particle changes 
when it crosses the boundary sites. The average density 
of particles is equal for the two fluids and has the same 
value for all simulations: 3.5 particles per site. We do not 
consider effects of wettability in the system, although it 
is possible to introduce it in the cellular automata model 
m- Thus, in the unstable situation studied below, the 
fingers are controlled by viscous forces at large scales and 
capillary effects at small scales. 

The bead-filling algorithm is as follows: first, we choose 
the total amount of sites to be solid sites, then we choose 
randomly (with an uniform distribution) the position of 
the bead center and its radius (between 20 and 50 lattice 
units). If the bead is sufficiently distant from any other 
solid site (20 lattice units is chosen), the new bead is ac¬ 
cepted and the total amount of solid sites increases. This 
process is iterated until the pre-defined total fraction of 
solid sites is reached or exceeded. This algorithm leads 
to the exponential distribution of radii shown in Fig. 
which is well fitted by the Poisson distribution and im¬ 
plies an average radius equal to 25 lattice units and an 
average porosity (the void fraction of the system) equal 



FIG. 7. (Golor online) Snapshots from three simulations with 
distinct parameters of the displacement of the blue fluid by 
the red fluid in a disordered medium, the scattering rate of 
the red (light gray) [blue] fluid is 0.008 [0.001] in (a), 0.001 
[0.001] in (b) and 0.001 [0.008] in (c). In (b) and (c) just below 
the image of the snapshot we present the interface between 
the two fluids. The size of the systems is 5000 x 2500 \/3/2 
and the forcing rate is 3.0 x 10“"*. In (b) the sites at the first 
and the last horizontal lines are solids. 
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to 0.81. 
(a) 




FIG. 8. (Color online) Height function h{y,t) for systems 
with size 10000 x 6000 \/3/2. In (a) the scattering rate of 
the red (light gray) and blue fluids are 0.001 (Case I), while 
in (b) the scattering rate of the red (light gray) fluid is 0.001 
and the scattering rate of the blue fluid is 0.002 (Case II). 
Time increases from left to right and each hgure represents 
one realization of the respective case. 

In Fig. [^we present snapshots of the displacement of 
a blue fluid by the red fluid in a disordered medium in 
three regimes, with the same forcing rate: 3.0 x 10“^. 
In Fig. the scattering rate of the red fluid is higher 
than the scattering rate of the blue fluid and the stable 
interface is flat in the void regions; in Fig. ^b) the scat¬ 
tering rate of the two fluids are equal and the interface 
remains compact; while in Fig. [^c) the defending fluid 
has a higher scattering rate compared to the invading, the 
interface displays a dendritic pattern, as expected from 
Saffman-Taylor instability, and shows bubble production. 

We shall now discuss in detail the interface dynamics 
in the following cases: 

• Case I: blue and red fluids have the same scatter¬ 
ing rate (0.001); 

• Case II: the scattering rate of the blue fluid is 
0.002, while the scattering rate of the red fluid is 
0 . 001 . 

A single-valued interface height function h{y, t) is defined 
as the horizontal position (with h increasing from left to 


right) of the rightmost red site in the horizontal line at y 
for the time step t. In Case I a constant forcing rate of 3 
X 10“"^ is applied, while in Case II a variable forcing rate 
ranging from 3 x 10“^ to 5.341 x 10“"^ is used, depending 
on the fraction of the system occupied by the red fluid. 
The height functions so produced are presented in Fig. 
for the two cases studied. 



t 



FIG. 9. (Color online) Interface width IT as a function of 
time for (a) Case I and (b) Case II for the indicated vertical 
dimensions, with a = '/Sl2, solid lines gives the best power- 
law function that best fits the last points of the curve for 
the largest system size. Insets: average value of the interface 
height as a function of time. 

The interface dynamics can be characterized by the 
time dependence of its width m- 

WiLy, t) = ^{[h{y,t)-h{t)f), (22) 

and the spatial dependence of the equal time correlation 
function: 

C(l,t) = yj{{h{y,t) - h{y',t)Y), for I = \y - y'l (23) 

where angular brackets denote averages on y and on dif¬ 
ferent realizations of disorder (our data were obtained 
from simulations of a minimum of 50 and a maximum of 
100 disorder realizations), while h{t) denotes the average 
value of h along y for a given disorder realization in time 
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step t. It is expected m that W and C presents, respec¬ 
tively, a power-law behavior for short time and length 
scales: 


W{Ly,t), for t tx; (24) 

0(1, t) for I (25) 

where /? is the growth exponent, a is the roughness ex¬ 
ponent, tx defines the crossover time above which 

the interface width saturates, where z is the dynamic ex¬ 
ponent, and ^ is the correlation length parallel to the 
interface. The crossover time tx is such that for t'^tx 
the interface width saturates in a value 

Wsat ^ Ly‘^- (26) 

It is worth mentioning that the ILG model presents an 
intrinsic fluctuating interface with a very short width 
and for system sizes with Ly < 64, the values a = 1/2, 
/I = 1/3, and z = a//3 = 3/2 were estimated [55], which 
are the exponents of the Kardar-Parisi-Zhang (KPZ) 
m continuum model. However, for system sizes with 
Ly ^ 64, logarithmic corrections were observed [55] in 
the scaling law for the saturation width. On the other 
hand, under a Saffman-Taylor instability, the interface 
width is not expected to saturate. 

In Fig. [^we present the growth of W and the average 
height < h > (see inset) as a function of time for the 
two cases studied. First, we notice that < h > is a linear 
function of time, as shown in the inset, and the velocity 
of the interface is constant, as mentioned above, with a 
value of 3.3 x 10“^ lattice units per time step in Case 
I and 3.6 x 10“^ lattice units per time step in Case II. 
The interface width exhibits two distinct behaviors with 
a crossover at t ~ 10000, and W « 50, which is approx¬ 
imately the average bead diameter: (i) for initial times, 
we find (3 « 0.53 and /3 « 0.51 for cases I and II, respec¬ 
tively; on the other hand, discarding the initial transient 
behavior, the exponent /3 is found to be /3 = 0.27±0.07 for 
Case I and /3 = 0.56 ± 0.07 for Case II. In Case I, for the 
large lattice sizes used {Ly = 1500-\/3/2 and 6000-\/3/2) 
and time interval of observation, i. e., (1 x 10^ - 2 x 10®) 
time steps, we have not witness any sign of width satu¬ 
ration, which may occur in Hele-Shaw experiments with 
quenched disorder [55] . 

In Fig. [^we present the equal-time correlation func¬ 
tions for the two cases at three values of time: ti = 
2000,t 2 = lOti and = lOOti. The curves present an 
upturn for short values of I, which is possibly related 
to the short scales fluctuations of the ILG model and a 
power-law behavior for intermediate values of I before 
saturation for I > 100; this length scale is higher for 
higher t and the value of W also increases with time. 
The difference in the density of scatterers in Case II im¬ 
plies a higher value for the the correlation at saturation 
and for the length scale associated with the crossover 
to this regime. We also calculate the time-dependent 
value of the roughness exponent a for 10 < / < 30, 



I 



I 



FIG. 10. (Color online) Equal time correlation function C{1, t) 
as a function of distance I for the indicated values of t in 
(a) Case I and (b) Case II. (c) Roughness exponent a as a 
function of time for the two cases. Simulations are performed 
in a system of size 10000 x 6000 \/3/2. 


which is shown in Fig. 10 (c). We see that a decreases 
slowly with time, starting from a « 0.5, and reaches the 
lower average value in Case II in the longer time studied: 
a = 0.39 ± 0.04 at t = 200000. 
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C. Universality classes of the roughening process 

In Case I, the estimated values of the exponents: /3 = 
0.27±0.07anda = (0.48 - 0.43)±0.03 for t > 2x10^, are 
very near the respective values for the discrete random 
deposition (RD) model with surface relaxation and 
the solution of the continuum Edwards-Wilkinson (EW) 
equation m, which are in the same universality class and 
have the following exponents [T^j: /3 = 0.25 and a = 0.5. 
The reasoning for this coincidence can be understood by 
noticing that in Case I the scattering rate of the blue 
and red fluids are equal, so that the interface dynamics 
is governed only by the random bead distribution and the 
surface tension. The effect of the beads distribution is to 
increase the width in an uncorrelated fashion (like the 
RD process), while surface tension softens the interface 
(surface relaxation effect), thereby inducing its correlated 
character. 



FIG. 11. (Color online) Same as in Fig. |8];b): Snapshot of the 
system as the most advanced portion of the red (light gray) 
fluid meets the right boundary of the lattice. 


In Case II the interface is unstable to viscous finger for¬ 
mation, in which case the viscosity contrast in Eq. (161 is 
A = 0.33. The simulation presents two important length 
scales, one of them is the average diameter of the beads, 
d « 50, and the other is associated with the viscous fin¬ 
ger dynamics, which is dominated by the perturbation 
with wavelength Am « 400, i. e., Ly/\m ~ 13. Com¬ 
paring Figs, [^a) and|^b), we notice a similar behavior 
for the height functions of these two simulations for early 
times; while, for later times, the height function shown 
in Fig. I^b) exhibits an oscillating pattern with wave¬ 
length ^ Am, superimposed by small spatial fluctuations 
with a characteristic length ^ d and caused by the effec¬ 
tive capillary pressure a jd [16) . As shown in Fig. 11 


for this time regime, the full integrity of the interface is 
somewhat lost and the system displays a huge produc¬ 
tion of bubbles, thereby originating a large mixed zone 
and an interface height function associated with the lim¬ 
its of this zone. In any event, the result of the quenched 
disorder manifests in a time dependent roughness expo¬ 
nent, shown in Fig. |10K c); in fact, by mapping the inter¬ 
face fluctuations at a fixed time onto a fractional random 
walk, with the y direction playing the role of the random- 


walker time, the roughness exponent is identified, after a 
short transient time, with the Hurst exponent of a subd- 
iffusive fractional Brownian motion |15j . i. e., a(t) < 0.5. 
More precisely, we find for Case I: a{t) = (0.55 - 0.43), 
while for Case II: a{t) = (0.53 - 0.39), with error bars in¬ 
dicated in Fig, [l^ c). Therefore, the self-afHne local frac¬ 
tal dimension |15L IM] . Df{t) = 2 — a{t), reads: Df{t) = 
(1.45 - 1.57) and Df {t) = (1.47 - 1.61) for cases I and 
II, respectively. In both cases, the early-time values of 
the fractal dimensions are consistent with the measured 
fractal dimension of immiscible displacement of viscosity- 
matched fluids in 2D porous media [SU], which was sug¬ 
gested [in] to be related to the fractal dimension of the 
external perimeter in invasion percolation [94) . In addi¬ 
tion, for Case II, it is quite satisfying to observe that our 
measured value of Df{t), at the breakthrough, is consis¬ 
tent with the fractal dimension of the percolation cluster 
backbone HQ], and also in agreement with fractal dimen¬ 
sions associated with: (i) fingering structure in a circular 
Hele-Shaw cell with quenched disorder and fluid invasion 
at the center of the cell for high Ca [66]; (ii) viscous fin¬ 
gering structure due to invasion of a less viscous fluid in 
a more viscous one in a horizontal Hele-Shaw cell with 
quenched disorder for side-sizes I (box counting) in the 
interval Ly > I > a/Ca [92] . 


IV. SUMMARY AND CONCLUSIONS 

In this work we presented a thorough study of a model 
based on lattice gas cellular automata techniques, includ¬ 
ing random dynamic scatterers. Our approach combines 
the best features put forward in this class of models: mi¬ 
croscopic kinetic rules, interface hydrodynamic fluctua¬ 
tions effects, and the presence of random dynamic scat¬ 
tering sites. To the best of our knowledge, the combi¬ 
nation of these techniques was not used in conjunction 
before the present work. In fact, the main achievement 
of our work is to show that the dynamics of the proposed 
model is in very good agreement with the observed dy¬ 
namics of immiscible two-fluid flows in Hele-Shaw cells. 

In order to fulfill this goal, we first showed that for 
a single fluid under a constant injection of momentum, 
the macroscopic Darcy’s law equation is recovered, and 
the ratio k/y can be finely tuned through the model pa¬ 
rameters. We then reported simulations of two distinct 
but related phenomena, (i) For two immiscible fluids 
with the less viscous fluid displacing the more viscous 
one, the model exhibits the Saffman-Taylor instability, 
in which case the number of fingers, shape, and time 
evolution is in very good agreement with the predictions 
from a linear stability analysis and experimental obser¬ 
vation of the interface dynamics. We also mention that 
the finger dynamics is quite complex; indeed, as time in¬ 
creases, some fingers can give rise to both bubble forma¬ 
tion and finger amalgamation, (ii) Next, we studied the 
dynamics of the interface between immiscible fluids under 
quenched disorder, simulated by randomly distributed 
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beads with random diameters obeying the Poisson dis¬ 
tribution. In the case of matched-viscosity fluids, the 
growth and roughness exponents present values suggest¬ 
ing the same universality class of the random deposition 
model with surface relaxation: /3 = 0.25 and a = 0.5. 
In the case of the less viscous fluid displacing the more 
viscous one, we find /3 ~ 0.5, and that the interface 
early-time dynamics is very similar to the former case; 
however, as time increases, viscous hngers develop with 
the subsequent production of bubbles in the context of 
a complex dynamics, thereby giving rise to an increas¬ 
ing mixed zone. Finally, we emphasize that the time- 
dependent behavior of the roughness exponents are quite 
similar in both cases, with a monotonic decreasing from 
their starting value: a ~ 0.5; indeed, after a short tran¬ 
sient time, this exponent can be identified with the Hurst 
exponent of a subdiffusive fractional Brownian motion. 
In particular, after the transient behavior, the early-time 
fractal dimension of the interfaces are consistent with the 
measured fractal dimension of immiscible displacement of 
viscosity-matched fluids in porous media, and with the 
fractal dimension of the external perimeter in invasion 
percolation. Further, in the case of a less viscous fluid 


displacing a more viscous one, our measured value of the 
fractal dimension at the breakthrough is consistent with 
the fractal dimension of the percolation cluster backbone, 
and in agreement with experimental results in Hele-Shaw 
cells with quenched disorder. 

The above summary clearly supports the conclusion 
that the two above-mentioned phenomena were success¬ 
fully simulated by our model, with insights and interest¬ 
ing features unveiled. We thus expect that the reported 
results will stimulate further theoretical and experimen¬ 
tal studies of two-fluid flow in porous media. 
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